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ABSTRACT 

The edges of magnetically-dead zones in protostellar disks have been proposed as 
locations where density bumps may arise, trapping planetesimals and helping form plan- 
ets. Magneto-rotational turbulence in magnetically-active zones provides both accretion 
of gas on the star and transport of mass to the dead zone. We investigate the loca- 
tion of the magnetically-active regions in a protostellar disk around a solar-type star, 
varying the disk temperature, surface density profile, and dust-to-gas ratio. We also 
consider stellar masses between 0.4 and 2 Mq, with corresponding adjustments in the 
disk mass and temperature. The dead zone's size and shape are found using the Elsasser 
number criterion with conductivities including the contributions from ions, electrons, 
and charged fractal dust aggregates. The charged species' abundances are found using 
the approach proposed by Okuzumi. The dead zone is in most cases defined by the 
ambipolar diffusion. In our maps, the dead zone takes a variety of shapes, including a 
fish-tail pointing away from the star and islands located on and off the midplane. The 
corresponding accretion rates vary with radius, indicating locations where the surface 
density will increase over time, and others where it will decrease. We show that den- 
sity bumps do not readily grow near the dead zone's outer edge, independently of the 
disk parameters and the dust properties. Instead, the accretion rate peaks at the radius 
where the gas-phase metals freeze out. This could lead to clearing a valley in the surface 
density, and to a trap for pebbles located just outside the metal freeze-out line. 

1. Introduction 

Trapping solid material is a basic requirement for forming planets enriched in the heavy el- 
ements with respect to the central star. One way to trap solid particles is through gas pressure 
gradients. If the pressure locally increases with distance from the star, radial force balance dictates 
the gas orbit faster than Keplerian. Embedded particles then feel a tailwind, receive orbital angular 
momentum from the gas, and drift away from the star, collecting near the point where the general 
outward pressure decline resumes. Such local radial pressure maxima can occur in various ways. 
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They can form at th e distance where the disk's m a in hea t source switches from accretion power to 



(jjohansen et al.ll2009l : lUribe et al 



2011 



stehar ihumination (IHasegawa &: Pudrital2010bllal. 120111). o r can be part of long-Uved zonal flows 



Pinilla et al 



20121 ). Pressure bumps can also be made by 



unsteadiness in the overall inward spiraling of the disk gas that feeds the growth of the central star. 
Material will pile up where the inflow becomes slower than average. 

The best cand idate to drive the disk's inflow is turbulence driven by the magneto -rotational 



instability or MRI (jBalbus fc Hawlevlll99ll : iHawlev et al.lll995l : iBalbus Hawlevlll998l ). Since the 



MRI relies on magnetic forces, the disk gas must be ionized enough to couple to the fields. The 
low temperatures rule out thermal ionization, but the disk surface layers are nevertheless partly 



( Gammie 


199fi: 


Glasseold et al 


1997) 


in the disk interior ( Sano et al. 


2000) 



2000), so that the grains' abundance regulates the flow of material 



through the disk. 

Time-unsteady flow leading to local pressure maxima can occur where the disk crosses the 
threshold from rnagnetically-coupled to de coupled, at the dead zones' inner and outer radial edges 
(jLvra et al.ll2009l : iDzvurkevich et al.ll2010l ). Away from the edges, pressure bumps could arise from 
changes in the rate of flow through the magnetically-active surface layers above and below the dead 
zone. Such interruptions in the la yered flow might come from a iump in the dust recombination 
cross-section across the snow line (iKretke &: Linll2007l : iBrauer et al.ll2008l ). 



The dead zone itself can participate in the overall radial flow. Even where the resistivity 
is high enough to damp away the MRI, and magnetic forces cannot directly dri ve turbulence, the 
stres s es can be non-zero due to waves propagating from n earby turbulent regions (jFleming &: Stone 



20031 : IDzvurkevich et al 



2011 



201C 



gas-p article interaction (IWeidenschilling &: Cuzzi 



1993 



Okuzumi &: Hirose l201ll) and to hydrodynamical i nstabilities and 



Johansen fc YoudinI 120071 : iJohansen et al 



). The dead zone stresses affect the timescales for material to pile up in orbit around the 
star, and so must be considered when attempting to use the ti me-evolving flow in layered disks to 



explain episodic accretion events such as FU Orionis outbursts (lArmitage et al 



2012 



Martin fc Livioll2012l : IZhu et al.ll2010l ). 



2001 



Martin et al. 



The size and shape of the dead zone thus govern many aspects of protostellar disks' evolution. 
The dead zone's edges are determined by the non-ideal terms in the induction equation corre- 
sponding to Ohmic, Hall and ambipolar diffusion. The Ohmic term is typically strongest in the 
dense disk interior, the Hall term at inter mediate densities and the ambipolar term in the disk's 



low-density atmosphere and outer reaches (jWardle fc Ng 



Perez-Becker fc Chiang 



2011 



1999 



Wardle 



20071 : 



Sano h Stone 



2002a 



The dead zone determ ined by the Ohmic diffusivity has been mapped under ioniz ation by cos- 
mic r ays (jSano et al. I I2OO0I ). X-rays (|llgner k Nelson|[2006l ) and stellar energetic protons (jTurner &: Drakd 
2OO9I ). Recently its extent was estimated with a more prec ise magnetic field distribution using a 
sequence of local shearing-box chemical-MHD calculations (IFlaig et al.ll2012l ). The impact of the 
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Ohmic diffusion might be reduced by the associated magnetic dissipa tion heating, at least n ear the 
star where temperatures approach the thermal ionization threshold (jMuranushi et al.ll2012l ). 

The Hall term has at most a weak effect on the saturated strength of the MRI turbulence, at 
least when t he Hall and Ohmic te rr ns are not too different in m agnitude, according to numerical 
calculations (ISano &: Stond l2002bl ) . IWardle &: SalmeronI (120121 ) however suggest that the results 
would differ at very large ratios of Hall to the other terms. When the Hall term is taken into 
account in this regime, the gas flow through the magnetically-active layers could be faster. 

Ambipolar diffusion becomes important when collisions are rare enough that the neutral gas 
drifts with respect to the plasma. An alysis of the MRI in the linear regime (jBlaes Sz Balbus 



19941 : iKunz fc BalbusI 12004 : 



DeschI l2004l ) shows that the instability is suppressed when the plasma 



transfers its momentum to the neu trals slower than the in stability grows. The non-linear evolution 
of the MRI was first mo deled by iMac Low et al.l (|l995l ) in two spatial dimensions, followed by 



Brandenburg et al.l (|l995l ) for the strong coupling regime. In 3-D two-fluid loc a l-box simulations 



the non-linear MRI under ambipolar diffu s ion \y as studied by iHawlev &: Stond (|l998l ) for several 
different magnetic field configurations. iBail (|201ll ) show that the turbulent accretion stress declines 
smoothly with increasing ambipolar diffusivity, from around 1% of the gas pressure at the threshold 
for MRI turbulence down to 7 x 10~^ at diffusivitie s ten times greater. The local shearing-box 
simulations were unstratified. In ISimon et al.l ()2012l ). MRI turbulence on azimuthal background 
magnetic fields was modeled for conditions similar to those thought to occur in the outer reaches 
of protostellar disks. The turbulence died out almost completely below the threshold ambipolar 
diffusivity. 

The ionization degree of the gas is affected by the dust surface area available to sweep up 
the ions and elec trons. The recombin a ,tion cross-section is o f ten treated as residing on compact 



spher ical grains (ISemenov et al.l l2004l : lllgner &: Nelson 



2006; 



Turner &: Sano 



20081 : 



Turner et al 



2OIOI') so that t he area decreases as the grains grow. Here we take the alternative approach of 



Okuzumi (j2009l ). supposing that at least the smaller grains bearing most of the cross-section remain 



as loose, fiuffy aggregates of smaller particles. The fractal dimension is close to two and so the 
area remains unchanged as the aggregates grow. Larger particles may be compacted by collisions 
or internal forces as the solid material evolves toward planet formation, but we assume the small 
fluffy aggregates continue to dominate the recombination cross-section. Like Okuzumi we neglect 
the electrical polarization that affects small compact grains' interaction with approaching charges, 
and compute the abundances of the ions and electrons, and the charge state distribution of the 
aggregates, from his governing equation (his eq. 30). 

The situation with the magnetic diffusivity is complicated by the fact that the grains carry 
charge. Small PAH-like particles, due to their low mass and re latively hi g h cyclotron frequencies, 
can signiflcantly change the disk's overall magnetic coupling (jBai 



2011 



2011 



Perez-Becker fc Chiang 



(Mohanty, Ercolano &: Turner, in press). More-massive compact grains have little direct 
effect on the electrical current because of their low cyclotron frequencies, even when charged to the 
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limit where further electrons are repelled by the electrostatic force lWardld (|2007l ) . Fluffy aggregates 
however, due to their large surface area, can accommodate many more elementary charges. We 
therefore include the aggregates' contributions to the conductiyities. 

In this paper we study the factors determining the dead zone's extent in the presence of fluffy 
grains. After presenting a fiducial model we vary the size of the monomers making up the fluffy 
aggregates, the aggregates' abundance relative to the gas, the overall disk mass and the surface 
density power-law index, the gas temperature, the stellar mass and the strength of the cosmic ray 
ionization. Across this parameter space we find the magnetically-dead region has a diversity of 
shapes and sizes. 

The paper is laid out as follows. In section 2 we describe how we calculate the ionization state. 
Section 3 has an outline of the magnetic dissipation mechanisms and the criteria for coupling the 
gas to the magnetic field. In section 4 we present the resulting dead zone shape as a function of 
the monomer size, dust-to-gas ratio and other parameters. We discuss the implications for building 
local pressure maxima in section 5, and summarize our findings in section 6. 



Fiducial Disk Model 



2.1. Gas and magnetic field 



We begin by mapping the dead zone in a fiducial disk model, taking stellar mass = IMq 
and normalizing the surface density to S = 1700 gcm"^ and the temperature to T = 280 K, both 
at 1 AU. Rather than a power-law surface den sity profile, we use a p ower-law with exponential 
outer roll-off, the similarity solution obtained by lHartmann et al.l ( 19981 ). This profile yields a good 
match between the radii of the disk measured in CO emission and in dust continuum iHuehes et al.l 



(|2008l ) . The surface density falls off with the radius r as 



exp 



where Sc, the surface density at the cut-off radius Tc, is equal to 

Mdisk 



(2-p). 



27rr2 



(1) 



(2) 



We take the median power-law index p = 0.9 observed by [Andrews et al.l (l2009l ) in the Ophiuchus 
star-forming region, and set the cut-off radius Vc to 40 AU. The surface density then matches the 
standard minimum-mass Solar nebula at 1 AU and again at 100 AU. The disk's whole mass of 
Mfjisk = 0.064Mq is stable against self-gravity, as the Toomre parameter is nowhere lower than 4. 
The density is easily reconstructed as 



27rcor sin(G) 



• exp 



cos 
"2cgsinG2 



(3) 
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where R = rsin(0) is the cyhndrical radius and cq = H/{rsinQ) the isothermal sound speed, 
with H the pressure scale height. Here we use spherical coordinates (r, Q, (p) with the rotation axis 
corresponding to G = 0. The orbital angular speed and temperature are 



^ = \T^^, and T = To./^i^, (4) 

with To = 280 K at Hn = 1 AU. The surface density profile with and without the exponential roll-off 
is shown in Fig. [TJ Cosmic rays reach the midplane largely unimpeded at r > 15 AU. Scattered 
X-rays similarly reach the midplane outside r > 60 AU. Nevertheless, the thin gas in the outer disk 
could be magneto-rotationally stable due to ambipolar diffusion. To map any stable regions, we 
need to obtain the densities of the charged species in a multi-component plasma. 

We calculate the conductivities in the domain extending from 1 to 200 AU in radius and over 
7r/2 ± 0.4 in the polar angle O. The disk's concave or flaring shape means the model covers ±12 
pressure scale heights at 1 AU but just itSH at 200 AU. 

We give the disk a mean vertical magnetic field, 

B = ^87rpc2//3p,o, (5) 

where the midplane plasma beta /3p^o = 400 is a constant and p{r) and Cs(r) are the gas density and 
sound speed at the midplane. The RMS vertical magnetic field in the saturated MRI turbulence 
varies only slightly with height up to the location w here the plasma beta falls below unity, in 



stratified local ideal-MHD shearing-box calculations bv lMiller &: Stond (|2000l ). To simplify matters 



we therefore make the field strength independent of height. 



2.2. Dust 



We assume a two-component dust population consisting of fractal porous aggregates and larger, 
more compact rocks. The fractal aggregates are assembled of small, nonaggregated du st (monomers) 
with radius an. The fluffy aggr egates couple well to the gas through drag forces (jMeakinl Il99ll : 
Blumll2004l : iDominik et al.l 120071 ) an d we assume they r emain well- mixed in the ga s, with dust-to- 
gas mass ratio fdn- Bot h numerical (IKempf et al.lll999l ) and experimental studies (jWurm &: Blum 



1998 



Blum et al.i ll998l : blum fc WurmI boool l show that such aggregates form at early growth 



stages when grain-grain collisions are slow. The outcome is an ensemble of fractal aggregates with 
fractal dimension D < 2, with a narrow, quasi-monodisperse mass distribution. The projected 
cross-section of the aggregate is u ~ (TqA^, where uo = ttuq is the geometrical cross-section of the 



monomer fsee lOkuzumil (120091') for det ails). The fractal growth continues until the aggregates reach 
centimeter sizes (jSuvama et al.ll2008l ). The compact bodies in contrast have D ~ 3. They may 
form in low-probability collisions between the fluffy particles, with impact speeds high enough to 
modify the aggregates' structure but too low for complete disruption. 



- 6 - 




Fig. 1. — Surface density and temperature profiles in the fiducial model disk. The dashed line shows 
the surface density power law with exponential roll-off. The solid black line shows the power-law 
alone, while dotted gray lines indicate the penetration depths of cosmic rays and X-rays, 96 gcm~^ 
and 8 gcm~^, respectively. The red line is the temperature profile. 
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The fluffy aggregates dominate tlie total cross-section for ions and electrons because of their 
much greater surface area per unit solid mass. The compact particles simply serve as a sink of 
solid material. The fluffy grains therefore determine the ionization state, even though the compact 
particles may contain most of the mass. 



2.3. Ionization and Recombination 



Thermal ionization is effective only within less than 1 AU of the central star, where tem- 
peratures above 1 OOP K mean alkali metals can ionize when colliding with hydrogen or helium 
(jUmebavashil 1 19831 ). At the lower temperatures found further out, ionization comes from cosmic 
rays, stellar X-rays and ultraviolet photons, and radionuclides. Our model includes all these except 
the ultraviolet photons, which penetrate a mass column less than 0.1 g cm~^ (Perez-Becker & 
Chiang 2011) so we do not consider them further. 



Most ionizations are of the abundant species, H2 and helium. In the total ionization rate 



C — CxR + CcR + C 



RN 



(6) 



the cosmic ray contribution is 

CcR = CcR^ 



exp 



D. 



3/4" 



DcrJ \DcrJ 



-4/3 



-4/3' 



(7) 



where the cosmic ray penetration depth Dqr = 96 gcm~ ^, and the mass column vertica lly above and 
below the point of interest are Da and -D^ respectively ([Umebavashi &: Nakandl2009l). The cosmic 
ray ionization rate in the in terstellar medium varies from 10~^® s~^ to 10~^^ s~^ (jUmebavashi &: Nakano 
198ll : iMcCall et al.ll2003al ). As a fiducial value we choose CcRO = 5 x 10 



-18 S-1. 



The X-ray contribution to the ionization rate is 



CxR 



Lx 



Xo29 erg s" 



Ci 

<2 



1 AU 
R 

+ e 



2.2 



(8) 



-15 



where the direct X-rays' intensity = 6x 10"^^ s~"^, the scattered X-rays' intensity C2 = 10~ s 
and the corresponding penetration columns A^i = 1.5 x 10^^ cm~^ or 0.006 gcm~^ and = 
7 X 10^^ cm~^ or 3 gcm~^ are measured by integrating the gas density vertically. This expression 
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for the X-r ay ionization rate was fit bv lBai &: GoodmanI (|2009l ) to Monte Carlo radiative transfer 
results from llgea &: Glassgoldl (j 19991 ). We take the young star's X-ray luminosity to be 10^^ erg s~^. 



Finally, we take a radionuclide ionization rate 

Cij7v = 7x 10-19 (/dg/10-2) 



(9) 



with the constant d erived from Solar system mea surements of the daughter products of the short- 
lived nuclide ^^Al (jUmebavashi fc Nakand l2009l ). Since the ^^Al is refractory and concentrated 
in the solids, and since the only solids well-mixed in the gas are the fluffy dust aggregates, the 
ionization rate is proportional to the aggregates' mass fraction f^g. 

Eqs. [MHl define the ionization rate throughout the cold parts of the disk. To find the ion- 
ization state we need a model for the recombination chemistry. A comple x network of chemica l 
reactions is involved both in the gas phase and on the surfaces of dust grains (jSemenov et al.ll2004l ). 
However in the presen ce of grains, the results can be well-approximated with a simplifie d network 
(jllgner &: NelsonI l2006l ) . We apply the simplified model put forward by lOkuzumil (|2009l ) . His key 
ansatz is that the grains have an approximately Gaussian charge state distribution. This holds if 
the grains are fluffy aggregates and (1) capable of holding large numbers of elementary charges, 
and (2) only weakly electrically polarized by approaching electrons and ions. The grain charge 
state distribution is then a function of a single dimensionless parameter F, the electrostatic energy 
between a charged grain and an incident electron, relative to the thermal energy. 

Writing the equilibrium rate equations for the electrons and ions, with the latter averaged 
over all ion species, and using electron-grain and ion-grain collision cross-sections averaged over 
the Gaussian grain charge state distribution, one can solve simultaneously for the ion and electron 
densities in terms of the same parameter F. Requiring charge neutrality then yields a single equation 
whose root is F. We solve this equation, Okuzumi's eq. 30, to obtain the number densities of the 
ions and electrons and the mean and width of the grain charge distribution. Our fractal aggregates 
have dust-to-gas mass ratio fdg and are made up of individual grains or monomers whose radius 
we choose between 0.005 and 5 /um. 

The most abundant molecular ion is often HCO+, since energetic particles convert H2 to H^ 
which quickly reacts with CO. However at temperatures high enough for metal atoms to desorb from 
the grains and enter the gas phase, such molecular ions readily transfer their charge to the metals. 
The metal ions recombine more sl owly and hence end up more abundant. As a representative metal 
we consider magnesium following lllgner fc NelsonI (j2006l ). 



The gas quickly reaches a charge equilibrium in the presence of the fluffy aggregates, with 
electron and ion densities 

C{md/mn) {VTT2^ - 1) ^^Q^ 



and 



fdgSeUeO- 5exp(-F) 



C{md/mn) {^/T+2^-l) 

fdgSiUiG g{l + F) 



(11) 
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where 

^ _ 2QCngas exp(r) 

Also Q is the mean gas-phase recombination rate, and the dimensionless parameter F = —{Z)e'^/{akBT) 
depends on the mean grain charge {Z) in units of the electron charge e, the grain radius a and ge- 
ometric cross-section a, and the temperature T. Other quantities are the electron and ion thermal 
speeds Ue,Ui and sticking probabilities on grains Se,Si for which we use Okuzumi's expressions. 

The grain charge state distribution is the Gaussian 



2(AZ2) 



(13) 



with spread 

1 + TakBT 

2 + r e2 

dOkuzumjboogl ). where the total dust number density ndust = pfdg /i^dust- Approximating the grain 



, = (14) 



charge state distribution by the Gaussian yields an overall charge-neutral plasma to high accuracy 
when the fluffy aggregates are made of at least a few hundred monomers each. In Appendix B we 
show the net charge that creeps in when the grains' charge spread is too small to be represented 
with a Gaussian. 

We compute the conductivity including the contributions from grains across the charge state 
distribution, along with the contributions from the electrons and ions. The free parameters are 
the individual grains' size and mass, the dust-to-gas ratio, the ionization rate and the gas-phase 
recombination rate. Our choices for the parameters governing the ionization state are listed in 
Table [TJ The gas-phase recombination is further discussed in the next section. 



2.4. Gas-Phase Recombination Rate 

A key weakness of this ionization recipe when applied over a large range of disk locations is 
that it does not capture the transition from one dominant ion to the next, and in particular between 
molecular ions and metals which occurs at the metals' thermal desorption temperature. The recipe 
yields a total ion abundance using a recombination rate averaged over all ion species x: 

rix, Ct = . (15) 

Here we focus on the transition between our representative metal ion Mg"*", with number density 
ni, and molecular ion HCO^, with number density n2. Under charge equilibrium, the metal ion 
number density is given by 

dni 

-Q^ = n2n^a - ciniUe - SiUiaudustni = 0, (16) 
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where is the number density of metal atoms remaining neutral and n[Mg] = ni + 71,3 the total 
number density of metal atoms in the gas phase. Also ci = 3 • 10"^^/\/T cm^ s~^ is the metal ions' 
radiative recombination rate, C2 = 3 • 10~^/\/T cm^ s~^ the molecular ions' dissociative recombi- 
nation rate and c.s = 3 • 10~^ cm^ s~^ the molecular-to- metal ion charge transfer rate coefficient 



(lllgner fc Nelsonll2006l ). Our task is to determine the effective recombination rate coefficient (eq. 9 
Okuzumil hoO^h ) by finding the ratio of metal and molecular ion densities y = ni /n2 in 



ct = {c2 + yci) . (17) 
1 + y 

First we determine the number density of thermally-desorbed metal atoms in the gas phase n[Mg], 
which depends on temperature. As described in Appendix A we balance adsorption of Mg on grains 
with thermal desorption. The metal-to-molecular ion ratio and gas ionization state are obtained 
by applying Okuzumi's recipe first assuming all the gas-phase metal atoms are ionized, and then 
iteratively correcting q to get the final y. 

In the fiducial model disk, the ion Mg^ proves to be the main source of free electrons at the 
midplane inside 5 AU. We introduce the concept of the metal line, the location at which a metal 
such as magnesium freezes out. In our fiducial model, the metal line lies entirely inside 8 AU where 
r~ 100 K (Fig. [21 top). 

In Fig. [2] (bottom) we show contours of the dimensionless factor g(T) introduced in eq. [T2j For 
^(r) ^ 1 gas-phase recombination dominates. The isolines of g{T) follow the gas density isolines 
in the outer disk (r > 50 AU), but in the inner disk this is no longer the case. Still it is safe to 
say that for gas densities p < 3 x 10~^^ gcm~^ the recombination occurs mostly in the gas phase 
rather than on the grains. In the outer disk, or in the inner disk above 5 pressure scale heights, 
the gas ionization degree is therefore insensitive to the dust properties. 



3. Magneto-Rotational Turbulence 



The gas couples to the magnetic fields, and gas movements can bend the field lines, if the fields 
diffuse slower than they are advected. For protostellar disks, various bending length scales and 
advection time scales have been suggested. If we take the length scale to be the gas scale height 
H, and the time scale to be the differential rotation or shear time 0"^, then the induction and 
diffusion terms in the field evolution equation, V x (i? x B ) and V xJrjW x B), are roughly QB 
and riB/H"^, leading to the condition cl/{-r]Vt) > 1 (eq. 35 in IWardlel (120071 )). 



However the turbulence involves structures smaller than the gas scale height. The magnetic 
field is sub-equipartition with the gas, so the Alfven speed va is less than the sound speed Cg 
and the MRI wavelength vaz/^ is less than the gas scale height Cs/^. Furthermore the poloidal 
field compo nents Bp, and Br are weaker than the azimuthal component in 3-D MRI turbulence 
calculations (jMiller fc Stonelboool ) . allowing the poloidal components to be bent into sharper curves. 



-9.0 -8.5 -8.0 -7.5 -7.0 -6.5 -6.0 




Radius [AU] 

Fig. 2. — Left panel: logarithm of the total gas-phase recombination rate coefficient q in the 
fiducial model after eq. [T71 in the (r, B — 7r/2) plane. Dotted lines mark where the temperature is 
150 K (left) and 100 K (right). Right panel: logarithm of g{T); gas-phase recombination dominates 
where g{T) ^ 1 (contour lines drawn in blue). 
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This is borne out by the observation that the turbuler it correlation lengths a r e anisotropic, wit h 
eddies reaching size H only in the azimuthal direction (jGuan &: Gammiell201ll : iFlock et al.ll2012l ). 



MRI turbulence under stratified conditions is probably sustained by the cycle in which (a) 
is buoyant and its rise generat es Bp,, (h) from .Bf) the MRI produces i?T-, a nd (c) the shear makes 
azimuthal field again from B^. (iTout &: Pringldll992l : iHirose &: Turnerll201ll ). Step (b) in particular 
creates magnetic structures with a short vertical wavelength and may be the first step to shut down 
as the diffusivity is increased. 

The appropriate criterion for sustained MRI turbulence is therefore that the fields diffuse 
across the MRI vertical wavelength in less than the MRI growth time. This corresponds to the 
dimensionless Elsasser number 



A = 



> 1, 



(18) 



where the Alfven speed contains only the vertical component of the magnetic field. We use the 
Elsasser number criterion through the rest of this paper. Numerical calculations with an Ohmic 
diffusivity show that weak MRI turbulence operates down to 



r?0 



0.1, 



(19) 



with larger, slowe r eddies and lower turbulent stresses (ISano &: Miyama 



1999 



Sano et al.l 120041 : 



Simon et al.ll2009l ). A similar condition holds for the ambipolar diffusivity (jBaill201ll ). We calculate 
the Elsasser number A in eqs. 1181 and 1191 including the Ohmic and ambipolar magnetic diffusivities, 
I] — + VA- We neglect the Hall diffusivity in the MRI turbulence criterion because numerical 
results show it ha s weak effects when not more than two orders of magnitude greater than the 
Ohmic diffusivity (|Sano Stonell2002bl ). 



In Fig. Owe present schematically the dead zone and MRI-active region, separated by a, tran 



sitional zone where 0.1 < A < 1. In classical layered disk models (e.g. lArmitage et al.l (120011 )) the 
disk consists of active and dead zones only. Here we further divide the dead zone into the tran- 
sitional and more-thoroughly dead parts. The transitional part will prove important for trapping 
solids near the dead zone's edge. 

An important consequence of using the MRI wavelength rather than the density scale height 
is that the resulting Ohmic dead zone varies in size depending on the magnetic field strength. Am- 
bipolar diffusion in contrast yields an Elsasser number independent of the magnetic field strength, 
since both the diffusivity and the squared Alfven speed are proportional to the magnetic pressure. 

We impose one more conditi on for MRI turbulenc e: since the instability grows quickly only 
on sub-thermal magnetic fields (jKim &: Ostrikerl |2000| ). we require that the plasma beta (3 = 

87rPgas/S2 > 1. 
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Table 1: Ionization parameters in the fiducial disk model 



Value 


Units 


Crn = 7 - 10-i3(/d9/10-2) 


s-i 


(cm = 5 ■ 10-1* 


s-i 


OixR) = 6 ■ 10-12 


s-i 


C2(XR) = 1 ■ 10-15 


s-i 


N = 400 




a = 20ao 


cm 


ao = 10-5 


cm 


fdg = 10-" 




(T = Tra^ 


cm^ 


(TO = Trag 


cm^ 
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Fig. 3. — Schematic representation of the MRI- active region (blue), the transitional region labeled 
'T' (yellow) and the magnetically dead zone (orange), based on evaluating the Elsasser number in 
eq. 27. 
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3.1. Magnetic Diffusivity 



The electrons, ions and charged grains all contribute to the conductivity, while coupling to the 
magnetic field in varying degrees. A charged particle x has cyclotron frequency 



\qx\eB/{mxc), 



(20) 



where g^e is the charge and rrix the mass. The damping time of the charged particle's motion 
relative to the neutrals is 

iTi -X- m . I 

(21) 



ITln + 



1 



nin Un < (TV >xn 

where {crv)xn is the momentum-transfer rate coefficient for colliding with neutrals, averaged over 
the Maxwellian velocity distribution. We adopt (av)pn = 1 0~^^ ■\/l28kfiT / (97rmp) for the electrons, 
and {(Tv)in = 1.9 x 10"^ for the ions (Draine et al. 19831 ). The corresponding rate coefficient for 
the dust grains follows eq. 21 in Wardle &: Ng| ( 1999 ). where the collision speed is the neutrals' 
thermal speed since the grains' drift through the gas is subsonic. Putting together the cyclotron 
frequency and the collisional damping time, we obtain the parameter hx = Tx'^x that describes how 
well the charged particle couples to the magnetic field. For electrons, typically 6e ^ 1- The quality 
of coupling to the magnetic fields drops with the mass of the particle. The ions are more than 
4 orders of magnitude heavier than the electrons, and the grains are orders of magnitude heavier 
still. 



We treat the magnetic diffusivity as outlined bv I War did (|2007l ). computing the electric currents 
from the number densities of charged particles obtained with Okuzumi's recipe. The current density 
in the partially-ionized gas is 



J = aoE,\\ + X + apE^ 



(22) 



wher e uo, cth-, op are the Ohmic, Hall and Pedersen conductivities ()Cowling)ll976l : IWardle Ng 
1993), 

ec v-^ 

(To = -^l^nx\qx\Ox, (23) 



OH 



op 



ec nxqxh\ 

X 

ec V - nx\qx\hx 



The induction equation 



dB 

'dt 

then has diffusivities 



— = Vx(iixS)-Vx r]oV X B + r]H{V x B) x B + r?^(V x B) 
at 



Vo 



c^aH , c^crp 
VH = - — ^ and r]A = - — 5- - Vo, 



(24) 
(25) 

(26) 
(27) 



with (T_L = J {ajj + CTp). 
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3.2. Accretion Rate 



Given the spatial distribution of the Elsasser number through the model disk, we estimate the 
mass accretion rate M{R) by summing the effective accretion rates within the magnetically-active, 
transitional and dead layers. 



M = 37rCsiJ(Sacaac + StOt + Sdeadttdcad), 



(28) 



where subscripts * ="ac", "T" and "dead" indicate the three respective layers, are the layers' 
mass columns and a^, are the corresponding stress parameters. 

There is ample evidence that MRI turbulence yields an accretion stress wjiff) proportional to 
the magnetic press ure, that is, agr ^ WRrh/Ppna oc as long as /3 > 1 and the gas dominates 
the total pressure (jHawley et alJ Il995l : iBai fc Stond l201ll ) . In verticall y-stratified disks, the stress 



parameter a increases with height as the local plasma beta decreases (jMiller &: Stond 12000 ) . and 
the turbulent layers closest to the disk's magnetically-dominated coro na transport significant mass 



(jFleming fc Stone 



2003 



Turner et al. 



201C 



below unity, the stress drops off quickly (jMiller &: Stond l200d ). 



Okuzumi fc Hirosdl201ll ). Once the plasma beta falls 



Since our magnetic pressure is independent of height, we make the stress uniform in the 
magnetically-active layer. We let the stress fall smoothly with decreasing plasma beta below a floor 
value Pf. 

ao if > /3/ 

/ 9 \ 

exp 



Oin 



For a fair match to shearing-box results we choose the floor value /?/ 



(29) 



active stress parameter ao is inversely proportional to the density, ao = 4oop(,^ e) 
the midplane gas density. 



1.85. The magnetically- 
where Pmidir) is 



Our stress parameter in the transitional layer ax behaves in the same manner with the disk 
height but is an order of magnitude weaker, having ax = aac/10. In the dead zone, a non-zero 
Reynolds stress leads to Odcad = Qn.mid/lO^i where ao.mid is the magnetically-active value of the 
stress parameter at the midplane (jPzvurkevich et al.ll2010l ). 



We investigate two approaches to the accretion stresses in disk regions with poor magnetic cou- 
pling. In t he first we simply make the stress parameter a step function ac ross the Elsasser number 



threshold ([Martin et al.ll2012l : iMartin &: Lividl2012l : lArmitage et al.ll200ll ). However, shearing-box 



MHD parameter surveys show a smooth roll-off in the stress at Elsasser numbers below unity 
(jSano et al.ll2004l ). suggesting this simple approach might not capture all the accretion flow. 



We therefore also consider a second approach in which the stress parameter varies smoothly 
with the Elsasser number. We use the uppe r boun d on a(A) found in unstratified shearing-box 
ambipolar-MHD calculations bv iBai &: Stond (|201ll . their Fig. 4). Focusing on their results with 
net vertical magnetic flux corresponding to plasma beta 400, and adding a term for the non-zero 
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dead zone Reynolds stress, we obtain the approximate fitting formula 

" = ttdead + «0 ^^^ _ ^^y (-^Uj 

where 

logio/(A) = -4 + AV3^, (31) 

and k = logio(A + 1). This results in a smoothly- varying turbulent stress parameter. The two 
approaches yield identical accretion stresses for Elsasser number A = 10, but at other levels of 
magnetic coupling, the step-function under- or over-estimates the accretion rate: M^^'^p > j\^sniooth 
for 1 < A < 10, and M"*°p < M^mooth f^^. lo < A < 100. Finally we require that the stress 
parameter not exceed 0.95, at which point transonic turbulent speeds would lead to rapid shock 
dissipation. 



4. Results 

In this section we present the dead zone calculated for the fiducial model disk and then vary 
the parameters in turn. The parameter study is summarized in Table [21 

The fiducial model's dead zone is shown in Fig. [J] along with a version in which the monomers 
are enlarged to 1 /im, keeping the dust-to-gas ratio fixed. The Elsasser number divides the disk 
into dead and MRI-active regions. We further divide the dead zone into a thoroughly-dead part 
and the "transitional" zone, where the Elsasser number is less than one order of magnitude below 
the threshold for turbulence. Consider first the top left panel in Fig. HI From the midplane to a few 
scale heights above, the Elsasser number changes by several orders of magnitude. By contrast the 
radial gradients are gradual, with the transitional region extending in the midplane from 8 to 20 AU. 
The conditions for MRI turbulence are satisfied at the midplane between 20 and 128 AU. The dead 
zone defined by Ohmic dif fusion alone, usin g the Ohmic Elsasser number Aq = v\^/ {rjo^)-, has its 



customary almond shape (jSano et al.ll2000l ). indicated in Fig. Hlby dotted lines where Ao = 0.1, 1 
and 10. When the ambipolar diffusion is also taken into account, the dead zone is much larger and 
includes a transitional region whose outer edge is shaped like a fish's tail. This shape is sensitive 
to the disk parameters. The fish's tail fins are made by the ambipolas diffusion's dependence on 
the ion density. There is a gentle dip in the ion density at 2 scale heights. If the ion densities were 
overall just 20% lower, the ambipolar diffusion would decouple the whole interior. The fish's tail 
shape is stronger still if the aggregates are made up of bigger monomers (Fig. IH top right). The 
bottom plots in Fig. |3]show the importance of treating the transition from metal to molecular ion. 
Using Mg"*" throughout leads to dramatically overestimating the magnetic activity. On the other 
hand, HCO"'' alone accurately describes the dead zone's outer edge but yields Elsasser numbers an 
order of magnitude too low in the MRI-active layers nearer the star. The transition between the 
two ions occurs at radii 4 < r < 9 AU, where the temperatures are 150 K to 100 K (see Fig. 2). 
This transition leads to the jump in the accretion rates. This happens because the Elsasser number 
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is roughly one order of magnitude larger when the metal ions are present in the MRI-active layers. 
Thus, the model with only Mg"*" ion produces an accurate estimate for accretion rate for r < 5AU, 
and the 'HCO"*" '-model should be used for the disk outside r > 9AU (see Fig. [25] in Appendix [A|). 



4.1. Characterizing Ambipolar Diffusion 

Figure [5] shows that ambipolar diffusion is more effective than Ohmic diffusion in damping 
MRI in our fiducial model and over a wide range of dust abundances. We therefore take a careful 
look at the criterion for MRI turbulence in the ambipolar diffusion regime. 

Ambipolar diffusi on has been characteri zed by two dimensionless numbers, the ion Reynolds 



number Rg = ^iPj/ O. (jHawlev Sz Stondl 19981 ) and the ambipolar Elsasser number A^ = vj^^/ {tja^) 
(jBai (fe Stonell201lh . The Reynolds number accounts only for the ions, while the Elsasser number 
includes the contributions of all charged species. The question we ask here is, where is Rci good 
enough, and where must the other plasma components be considered? 

In Appendix [C] we show that the two numbers are equivalent if ions couple to the field and 
dominate the Pedersen conductivity, which is larger than the Hall conductivity. In Fig. [6] we 
therefore plot the ion coupling parameter bi, the ions' contribution as a fraction of the total Pedersen 
conductivity, the ratio of Hall to Pedersen conductivity, and the ratio of the electron number density 
to the ion number density, in the fiducial model disk. The ion Reynolds number is unreliable in 
two regions. The first is the dense interior, where collisions thoroughly decouple the ions from the 
fields. The second is the disk atmosphere above 6H. There, gas densities are so low that grains can 
settle at supersonic speeds. This would make the grains the biggest contribution to the Pedersen 
conductivity. However since the gas drag and density profile are determined by heating and cooling, 
and potentially also by magnetic support, none of which we model in detail, we choose to neglect 
the settling for all our figures. Fuller treatment may be worthwhile in future. 

Fig. [7] is a map of the ratio Rci / Aa in the fiducial model disk. The two differ by a factor 
more than 1000 near 1 AU in the midplane. However the ratio is near unity in the large region 
shown in light gray, which includes the boundaries of the MRI-active region. The black dashed 
lines marking where Rei = 1 almost overlap the yellow lines indicating where A a = 1- Yet the 
curve where Rei/ Aa = 1 is in quite a different location. This is because even within the light-gray 
region the two numbers differ slightly (less than 1%). 

The map of the MRI-active region looks basically identical whether drawn using the ion 
Reynolds or ambipolar Elsasser number. In our fiducial model, and in all models from Table O the 
ion Reynolds number is as good as the ambipolar Elsasser number for estimating the size of the 
MRI-active region. However, far from the active region's edge, the ion Reynolds number can be 
unreliable as an estimator for the turbulent stress. 
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Fig. 4. — Dead, transitional and MRI-active zones for the fiducial model (left column) and a 
version with larger monomers = 1 f-im (right column). Blue regions are MRI-active, yellow 
are transitional and orange are MRI-dead. The top row includes the transition between the Mg"*" 
and HCO"'' ions, the middle row has Mg+ everywhere and the bottom row has HCO"'' everywhere. 
Black solid contours indicate where the overall Elsasser number is 0.1, 1 and 10. Black dotted 
contours are the corresponding surfaces for the Elsasser number computed from Ohmic diffusivity 
alone. Black dashed lines mark heights H and 5H. 
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Fig. 5. — Contributions to the total Elsasser number (black dashed curves) from Ohmic diffusion 
(purple), Hall effect (green) and ambipolar diffusion (red) at radius 7.5 AU in the fiducial model disk 
and versions with dust-to-gas ratios from 10~^ (top left) to 10~^ (bottom right). The horizontal 
axis is the height in units of the pressure scale height, with cq = c^/uKepier- 
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Table 2: Sets of models in the parameter study. The rows are the monomer radius ao, the dust-to- 
gas mass ratio in the fluffy aggregates fdg, the disk mass in units of the stellar mass M^/M^,, the 
surface density power-law index p, the temperature at 1 AU Tq, the stellar mass in Solar masses 
M^:/Mq, the cosmic ray ionization rate outside the disk CcBO (eq- 6), and the ionization rate due 
to radionuclides (rn (eq. 9). 
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Fig. 6. — Top left: ion coupling parameter bi versus height at four radii in the fiducial model disk. 
Top right: ratio of ion Pedersen conductivity to total Pedersen conductivity. Bottom left: ratio of 
absolute Hall conductivity to Pedersen conductivity. Bottom right: ratio of electron to ion number 
densities. 
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Fig. 7. — Gray contours and gray dotted lines show the ratio Rei/ Aa in the fiducial model disk. 
The red line marks iJe^/A^ = 1. On the black dashed line Rci = 1 and on the nearly-overlapping 
yellow solid line Ka = 1, two estimates of the border of the MRl-active region. The green dotted 
lines are logarithmically-spaced isodensity contours labelled in gcm""^. 
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4.2. Magnetic diffusivities in the ion-dust limit 



Ambipolar is the largest of the three diffusivities near the midplane at dust-to-gas mass ratio 
10^^ (Fig. [5|). This is unusual, as Ohmic diffusion typic ally dorninate s near the midplane for 



weak magnetic fields and Hall diffusion for stronger fields (jWar did 120071 ). We therefore examine 
the reasons. The ionization balance falls in the ion-dust limit, where the grains sweep up many 
electrons and ions, becoming the most abundant plasma component with a mean negative charge 
of a fraction of an electron. Even so, the grains' large mass means they couple very poorly to the 
magnetic fields, leaving the Ohmic and Hall diffusivities dominated by the well-coupled electrons, 
and the Pedersen conductivity by the poorly-coupled ions: 



The Ohmic conductivity comes from the electrons because they couple b est to the field s; 
gq oc nh in eq. [251 and h^. ~ 10006j, while in the ion-dust limit ~ nj/100 (jOkuzumil |2009| ) . 



The Hall conductivity also comes from the electrons since ou oc n in eq.m The 



coupling 6e ^ 1 ^ &i means the electrons' contribution ~ Uf, is greater than the ions' ~ nj6? 



The Pedersen conductivity is in the ions because they are more numerous than the electrons, 
while the two coupling parameters bracket the maximum at h = 1 in cjp oc nb/(l + b'^) (eq. [ 



Consequently the Ohmic conductivity is greater than the Pedersen, which in turn is greater than 
the Hall, since {ao ~ Uebe) > {crp ~ riibi) > {an ~ nj/100). Thus in our eq. [271 o"± ~ and 
the ambipolar diffusivity rjA oc cjp^ is bigger than the Ohmic diffusivity rjo oc ctq^ . The special 
circumstance leading to ambipolar dominance in the ion-dust limit is our chosen magnetic field 
strength, which couples the ions too little for them to contribute to the Hall conductivity, but 
enough to dominate the Pedersen conductivity. 

For fdg < 10^^, the reduced grain surface area means the midplane near 7 AU no longer 
reaches the ion-dust limit. The ionization fraction is larger, the electron density is comparable to 
the ion density, and both exceed the density of the charged grains, which carry many electrons on 
average. The Ohmic and Hall conductivities are both larger in proportion to the increased electron 
density, while the Pedersen conductivity scales with the lesser increase in the ion density. As a 
result the Hall overtakes the Pedersen conductivity, giving the ordering ao > cth > crp- Since 
rjH ~ l/o"// > ~ ^p/^'h-, the largest diffusivity is the Hall one. 

The Hall term is special in that it depends on the sign of the magnetic field, suppressing the 
MRI for one field direction and amplifying it for the other. It is unclear whether the Hall term 
has a strong impact on a turbulent eddy penetrating into the transitional or dead zone, as an eddy 
typically carries magnetic field of both signs. The importance of the Hall effect in the nonlinear 
evolution of the MRI remains to be clarified with further numerical calculations. Therefore we here 
confine ourselves to considering the Ohmic and ambipolar Elsasser numbers. 
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4.3. Magnetic Field Strength 

One of the greatest uncertainties in modeling the magnetic coupling in protostellar disks is 
the strength and configuration of th e magnetic field. A lower limit on the field strength is the 



mean Galactic field of ~ 10 yuGauss (jBeck &: Hoerneall996l : iBeck et al.lll996l ). An upper limit is 



equipartition of the vertical magnetic field with the gas, which makes the shortest unstable MRI 
mode longer than the disk thickness. In our models the midplane plasma beta is independent of 
radius, while the magnetic pressure is independent of height (see eq. 5 in section 2.1). 

The ambipolar Elsasser number is basically independent of the magnetic field strength, but 
the Ohmic Elsasser number is proportional to (eqs. 21, 25 and 27). Fig. 5 shows that the 
ambipolar Elsasser number is smaller than the Ohmic Elsasser number only by a small factor near 
the midplane. Thus in a magnetic field just a few times weaker, the MRI would be damped by 
Ohmic rather than ambipolar difi'usion near the midplane. Another concern is above what height 
the disk is magnetically-dominated, with plasma beta less than unity. To clarify these points, we 
vary the midplane plasma beta from 10 to 10^. The corresponding magnetic fields are between 
3.25 and 10"^ miUi gauss at 100 AU, where th e field strength can potentially be measured using 



sub-millimeter polarization (jKreinv et al.ll2009l ). In Fig. [5] (left) we plot the radial locations where 
the Ohmic and ambipolar Elsasser numbers are 0.1 and 1, as functions of plasma beta. This range 
corresponds to the dead zone's outer edge, which we call the transitional zone. Its width shows how 
gradual the transition is between dead and active zones. We see that the Ohmic dead zone grows 
with plasma beta. At /3 > 4000, Ohmic rather than ambipolar dissipation sets the outer dead zone 
boundary in the midplane. On the other hand, for strong magnetic fields such as 3 milligauss, or 
midplane plasma beta 10, the Ohmic dead zone lies deep within the ambipolar dead zone. 

Choosing low midplane plasma beta also leaves the disk's upper parts strongly magnetized 
and stable against MRI. In our fiducial model, the plasma beta is below unity above 3.5H (Fig. [8] 
right, green line). For a field strength at 100 AU matching the galactic field of 10 /xG (plasma 
beta of 10^), the magnetized region is only above 5H (Fig. [HI right). Note that ambipolar diffusion 
defines the dead zone's minimum radial extent. A larger dead zone indicates that the fields are 
weak enough for the Ohmic diffusion to control the evolution of the MRI. 



4.4. Monomer Size 

Here we discuss the influence of monomer size on the shape of the MRI-active region (set 1, 
Table [2]). The outer border of the active zone is well aligned with the gas isodensity surface 
p = 10"^^ gcm~^. The active zone encloses an interior dead and transitional zone shaped like a 
fish tail. Outside 10 AU, the active zone is defined by ambipolar diffusion. The Ohmic diffusion 
plays hardly any role in the radial extent of the dead zone for the assumed magnetic field strength. 
However Ohmic diffusion is important for the dead zone's vertical thickness in the inner disk. The 
transitional layer widens outside 8 AU, giving the fish-tail shape. The widening is caused by the 
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Fig. 8. — Left: The radial extents of the Ohmic and ambipolar transitional layers at the midplane 
in the fiducial disk model, as functions of the plasma beta parameter. Gray shading indicates where 
the ambipolar Elsasser number is between 0.1 and 1. Diamonds show Ohmic Elsasser numbers of 
0.1 (bottom chain) and 1 (top chain). Right: Contours of plasma beta as a function of the height 
in pressure scale heights (horizontal axis) and the midplane plasma beta (vertical axis). On the 
right-hand vertical axis is the strength of the corresponding vertical magnetic field in Gauss at 
100 AU. A green line marks our chosen midplane plasma beta parameter (3 = 400. 
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switch from Mg^ to HCO"*" as the main ion, leading to lower ionization fractions since the molecular 
ion recombines faster. For very small monomers, 0.005 (not shown) and 0.01 /im in radius, the 
active layers are thinner and separated from a midplane island of magnetic activity in the outer 
disk (Fig. O bottom). Also, the ratio Rei/ Aa deviates strongly from unity. Small grains couple 
better to the magnetic fields and participate in currents in the disk corona. Meanwhile in the 
interior dead zone, the small grains readily sweep up free electrons, leaving rij ^ rig. These results 
suggest that the effects of very small grains such as PAHs are best characterized using the Elsasser 
number criterion, which includes the contributions of all the plasma components as outlined above 
in section 14. li 

The transitional zone gains in radial extent at the midplane as the monomer size decreases. 
For oq = 1/um, the transitional zone runs from 6 AU where A = 0.1 to 15 AU where A = 1, giving 

9 AU of transition between adead ~ 10~^ and Oac ~ 10~^ (section 3.2). For ao = 0.01/im, the 
transitional layer extends from 19 to 45 AU, a width of 26 AU. 

The accretion rates are plotted against the distance from the star for the various monomer 
sizes in Fig. [TO]using the approach set out in Section 3.2. Very small monomers with qq ^ 0.01 /im 
produce accretion rates below 1Q~^Mq yr~^ at 1 AU. The biggest monomers, ao = 5//m, give a 
maximum in the accretion rate between 7 and 8 AU. At about 10 AU the accretion rate has a 
gap, related to the radial variation in the active layers' thickness shown in Fig. [9l When the stress 
parameter a(A) is given a step-function shape, the accretion rate jumps from 9 x 10~^Mq yr~^ at 
7 AU down to 2 x 10~® at 12 AU. Aggregates composed of the intermediate-sized monomers yield 
jumps almost as large: from 6.5 to 2 x lO~^M0yr~^ for oq = 1 //m, and from 5 to 2 x 10"*^ for 
ao = 0.5 fim. 

However the accretion rate profile changes drastically when the smooth function is adopted 
for Q;(A). The accretion rate peaks near the metal line, and accretion is much slower in the outer 
disk. The reason is that the Elsasser number is very near unity along the whole midplane outside 

10 AU, so that the smooth function yields a turbulent stress that varies slowly across the border 
between transitional and active zone. 

4.5. Dust-to-Gas Ratio 

In this section we vary the dust-to-gas ratio (set 2, Table [2]), leaving other parameters un- 
changed from the fiducial model. The resulting dead zones and the ambipolar Reynolds-to-Elsasser 
number ratio are shown in Fig. [TTl The radial extent of the transitional layer again shows striking 
variations. The transitional layer's midplane radial width is 25 AU, 12 AU, 7.5AU and 2 AU with 
dust-to-gas ratios 10""^, 10~^, 10~^ and 10"^. From the shift in the dead zone shape with increasing 
dust depletion, we see that the fish-tail shape appears when the active surface layers merge with 
the midplane active zone (Fig. [TTl top left and top right). At the lowest dust-to-gas ratio, the 
transitional layer splits into two pieces, located near 1.5 and 10 AU, and both about 2 AU in radial 
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Fig. 9. — Contours of the total Elsasser number (black solid lines) and ohmic Elsasser number 
(black dotted lines) in disks containing fluffy aggregates composed of monomers of various radii. 
The MRI-active regions are shown in blue, the transitional regions in yellow and the dead zone in 
orange. 
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Fig. 10. — Accretion rate in solar masses per year for various monomer sizes. The dotted line shows 
the accretion rate for the case when the whole disk is active with the turbulent stress from eq. 29. 
Note that the smooth- function turbulent stress can lead to higher accretion rates in the A > 10 
regions (see section 3.2). 
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extent at the midplane (Fig. Illl bottom right). The MRI-active region's outer boundary again hes 
at the gas isodensity surface p = 10~^^ g cm~'^. 

The radial profiles of the accretion rate are shown in Fig. [T2j For the step function a (A) and 
dust-to-gas ratios less than 10~^, the accretion flow shows a gap from 10 to 12 AU. The deficit in 
the gap is as much as a factor of four for dust-to-gas ratio fdg = 10~^. However no such feature 
appears with the smooth a(A) function. Instead the accretion rate steps up inside the metal line, 
located between 5 and 9 AU, and the dead zone edge has no clear signature in the radial accretion 
rate profile. 



4.6. Disk Mass 

In set 3 from table [21 we vary the disk mass by scaling the fiducial model's densities from eq. 3 
(section 2.1). Increasing the mass makes more of the outer disk active, if we hold fixed the ratio of 
gas to magnetic pressure (Fig. [T3l) . The reason is that higher densities make the ambipolar diffusion 
term less important relative to the induction term. The outermost border of the active zone shifts 
from 100 AU for the disk with Mdisk = O.SMfiduciai to 140 AU for Mdisk = l-SMfiduciai, following 
the gas isodensity contour p = 10"^^ gcm~^. At the same time, increasing the disk mass has little 
effect on the location and thickness of the active layers at higher latitudes. The active layer base 
simply moves slightly further from the midplane, as does the layer absorbing the ionizing cosmic 
rays and X-rays. The active layer becomes a little thinner because its top is still constrained by 
the requirement that magnetic pressure is not too much greater than gas pressure. 

For the smallest disk mass in set 3 (0.2Maduciai) MRI-active layers are found only at the higher 
latitudes, and the midplane is either dead or transitional. For larger disk masses, the midplane 
active island appears (Fig. [13l top right). In the two most massive disks from set 3, the active 
surface layers and midplane island merge. Over the same mass range, the interior dead zone grows 
due to the reduced ionization at the higher mass columns. 

The transitional region between active and dead zones generally narrows with increasing disk 
mass, as the expanding inner dead zone approaches the expanding outer active island. The transi- 
tion region is 13 AU, 12 AU and 10 AU wide for Mdisk = 0.5, 0.8 and l.SAfgduciai- At the lowest disk 
mass, the transitional region stretches in the midplane from 4 to at least 200 AU. A substantial 
volume is thus capable of MRI turbulence at reduced levels. 

The radial profiles of the accretion rate are shown in Fig. [141 Note that we keep the dust-to-gas 
ratio and monomer size same as in the fiducial model. The accretion rate is about 10~^ yr~^ 
at our inner boundary. For the smooth a{K) function the accretion rate dips at 9-10 AU near the 
metal line. In the more massive disks, the accretion profile becomes less dependent of mass. 

The radial profiles of the accretion rate are shown in Fig. [TH Note that we keep the dust-to-gas 
ratio and monomer size the same as in the fiducial model. The accretion rate is about 10~® yr~^ 
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Fig. 11. — MRI-active regions appearance for various dust-to-gas ratio, see also Fig. [9j Colors and 
symbols are as in figure [9j 
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Fig. 12. — Accretion rate in solar mass units per year calculated for the MRI-active sur- 
face density and for /? > 1 disk regions, for disks with various dust-to-gas ratios, log{fdg) = 
[—2, —3, —4, —5, —6, —7]. The dotted line shows the accretion rate for the case when whole disk is 
active with the turbulent stress taken from eq. 29. Note that the smooth-function turbulent stress 
can lead to higher accretion rates in the A > 10 regions (see section 3.2). 
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Fig. 13. — MRI-active regions' appearance for disk mass 0.2, 0.5, 0.8 and 1.5 times the fiducial 
model. Colors and symbols are as in figured! 
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at our inner boundary. For the smooth a(A) function the accretion rate dips at 9 AU near the 
metal hne. In the more massive disks, the accretion profile is almost independent of the mass. 



4.7. Surface Density Slope 



In set 4 (Table [2]) we vary the power-law index p in the surface density profile E oc from 
0.5 to 1.5, keeping the disk mass and other para meters as in the fiducial r nodel. The slope p = 0.5 
is frequently adop ted for numerical simulations ([Fromang &: Nelsonll2006l : iDzvurkevich et al.ll2010l : 



Flock et al 



2011 



Milli meter interferometer measurements of the dust continuum emission yield 
a median index p = 0.9 (jAndrews et al.ll2009l ). while Hayashi's minimum-mass Solar nebula model 
has p = 1.5 to approximate the distribution of mass in the Solar system. 

With increasing surface density slope, the active zone expands outward from 100 AU to 195 AU 
(Fig. [T5]) . At the same time the midplane dead zone thickens near 1 AU. For slopes p between 0.5 
and 1.3 we obtain a single MRI-active region with a fish-tail shaped transitional layer. For p = 1.5 
the outer disk's low densities make the ambipolar diffusion stronger still. The transitional region 
is larger and isolates the magnetic activity in the two surface layers from the midplane island of 
activity. The transitional layer between the island and active layers at 50 AU (Fig. [T5l bottom 
right), is only about one s cale height thick. As MRI turbulence can overshoot the dead zone 
boundary by a scale height (jTurner et al.ll2010l ). the transitional layers could be turbulent thanks 
to the activity above and below. 

In Fig. [16] we plot the accretion rates for the models in set 4. The inner border of the MRI- 
active region moves inward very little as the slope is increased. Similar to the results in the previous 
section, the accretion flow shows no gaps under the step- function a{A). The accretion rate is about 
10-8 M0 yr-i at our 1 AU inner boundary. For the smooth a (A) function, M peaks near the metal 
line. The accretion rate dip is deeper for shallower surface density profiles. 



4.8. Temperature 



Here we explore Set 5 (Table [2|) in which we vary the gas temperatures. The domain we con- 
sider remains too cool for significant thermal ionization. Nevertheless the temperature is important 
for the conductivity because it regulates the disk thickness and hence the gas density, determines 
the gas-phase metal abundance, and sets the particles' thermal speeds. We consider a range of 
tempera tures from those observed at the starlight -heated surfaces of disks in nearby star-forming 
regions (jBouwman et al.l l200ll : iFurlan et al.l l2009l) down to those inferred for the disk midplanes 
through radiative transfer modeling ([Chiang &: GoldreichI [19971 : [Pullemond et al.l [20071 ) . We con- 
sider temperatures at 1 AU, Tq, between 150 and 600 K, keeping the temperature's radial power-law 
slope constant at —0.5. 
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Fig. 14. — Accretion rate vs. radius in set 3 where the disks have masses 0.2, 0.5, 0.8, 1, 1.2 and 1.5 
times the fiducial model. Dotted lines show the accretion rates if the disks were active throughout. 
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Fig. 15. — MRI-active regions appearance for surface density slope 0.5, 0.9, 1.3, 1.5, assuming a 
constant disk mass. Colors and symbols are as in figured! 
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Fig. 16. — Accretion rate vs. radius in the disks with surface density slopes p = 0.5,0.9, 1.3 and 
1.5. Dotted hnes show the accretion rates if the disks were active throughout. 
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The resulting dead zones are shown in Fig. [T71 In each case the MRI-active layers contract 
across the radius where the metal ions freeze out. The magnesium in our models freezes out at 
2.5 AU if To = 180 K, 4.5 AU if 280 K and 11 AU if 350 K. At Tq = 450 K the MRI-active layers 
reach deep towards the midplane. We observe that the transitional layer's fish-tail shape is present 
only for Tq < 450 K. For still higher temperatures, simple almond-shaped dead and surrounding 
transitional zones appear. The remnants of the fish-tail appear as islands in yellow roughly 2H 
above and below the midplane (Fig. 1171 bottom). Although the islands stretch over several tens of 
AU in radius, their vertical thickness is less than one scale height, and they could conceivably be 
activated by the neighboring better-coupled regions. 

The accretion rates for the models in set 5 are shown in Fig. [THJ The islands of transitional 
coupling for temperatures Tq from 450 to 600 K, seen in yellow in the bottom panels of Fig. \T7\ lead 
to small reductions in the accretion rate if we take a(A) to be a step function. If we instead use the 
smooth a(A) function, the accretion rate peak shifts outward following the isotherm T « 100 K, 
reaching 50 AU for the hottest disk (Fig. [18l bottom) . 



4.9. Stellar Mass 

In set 6 (Table [2]) we vary the stellar mass, which determines the orbital frequency, and through 
the vertical component of gravity affects the disk's gas density. Furthermore the stellar luminosity 
grows with the stellar mass, affecting disk temperatures. We let the disk mass be proportional to 
the stellar mass, Mjisk = 0.064M*. Luminosities for sta rs between 0.3 and 2 M(7) are taken from 



evolutionary tracks computed with the STELLAR code (jBodenheimer et al.l 120071 ). For simplicity 



we neglect the influence of the gas accretion on the evolutionary tracks. 

The stellar luminosity is plotted vs. the stellar mass in Fig. [T9j We choose the time-step closest 
to 1 Myr as a typical age for disk-bearing stars. From the luminosity we reconstruct temperatures in 
the circumstellar disk as follows. The location of the dust sublimation front, where the temperature 
is about Tg = 1500 K, can be calculated as follows: 

^•" = ^-(|)'=(4^)'"- <^^> 

From this point we assume the temperature falls off as T = Tgy^Rin/R) corresponding to the 
optically-thin approximation used in the minimum-mass Solar nebula model on which our fiducial 
model is based. This overestimates the temperature since the disk is optically-thick, but our point 
here is to see whether the dead zone depends on the temperature. 

For masses 0.4, 0.8, 1.2, 1.6 and 2 Mq, the temperatures at 1 AU are 446.8, 518.7, 559.3, 
604.6 and 650.9 K, respectively. The Solar-mass star aged 1 Myr has a higher luminosity than our 
fiducial model, yielding a temperature at 1 AU of 532.3 K instead of 280 K. 



Higher luminosities stretch the metal-dominated region further from the star than in the 
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Fig. 17. — Magnetic coupling maps for disks with a range of temperature profiles T 
To(r/lAC/)"°•^ with To = 180, . . . , 500 K. Colors and symbols are as in figure SI 
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Fig. 18. — Accretion rate vs. radius in disks with various temperature profiles T = T^ir jXAU^ 
where Tq = 180, . . . , 600 K. Dotted fines indicate tfie accretion rates when the whole disk is active. 
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Fig. 19. — Stellar luminosity vs. stellar mass from evolutionary tracks at 1 Myr. A star symbol 
indicates the stellar luminosity in our fiducial model. 
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fiducial model (Fig. [20|) . The Elsasser number A = 10 contour stretches along with the Mg"*"- 
dominated region, reaching 20 AU for the 0.4-Mq star and 45 AU for the 2-Mq star. When the 
metal-dominated region reaches much beyond 10 AU, the fish-tail shape no longer appears in the 
transitional layer. Disks around stars with masses 0.4 to 1 Mq (top panels in Fig. [20|) have once 
again islands of transitional coupling, located about 2H above and below the midplane (see sec- 
tion 4.7). The threshold between Mg"*" and HCO"'' makes itself visible in the shape of the MRI-active 
region's outer border too. For the disks around stars from 1.2 to 2 Mq, one can clearly see the 
MRI-active Mg"*" region overlapping with the HCO"*" region (bottom panels in Fig. I20|) . Interest- 
ingly, while the Mg"*" -dominated layers extend further out around heavier stars (see for example 
the A = 10 contours), the ambipolar term still beats out the induction term at r > 40AU and 
3 > (0 — 7r/2)/co > 2. As this effect occurs outside the exponential surface density cut-off, where 
the density falls rapidly with radius, we expect this would have little effect on the accretion rate. 

The disks around the lower- mass stars have larger aspect ratios H/R or c^/uxep- The Keplerian 
velocity UKep increases faster with the stellar mass than does the sound speed. Since we fix the 
ratio of disk to star mass, the surface density increases with the mass of the central object. This 
leads to larger dead zones in the more massive disks. We observe in Fig. [20] that indeed the disks 
around the higher-mass stars have smaller MRI-active regions. The separation between the inner 
and outer radial boundaries of the transitional layer on the midplane changes slowly with the stellar 
mass, being [6, 9, 10, 11, 12, 13] AU around stars of [0.4, 0.8, 1, 1.2, 1.6, 2]Mq . Thus over the stellar 
mass range, the transitional layer widens from 6 to 10 AU, while the transitional layer's inner edge 
{log{h) = — 1 contour, or yellow-to-orange transition in Fig. [20]) moves outward by 4 AU. For each 
stellar mass the accretion rate peaks near the Mg"*" freeze-out radius, beyond which the outer disk is 
less strongly accreting, considering the smooth dependence of the stress parameter on the Elsasser 
number (Fig. 121] bottom). 



4.10. Cosmic Rays 



The final parameter we vary is the cosmic ray ionization rate. We consider values 10, 100 
and 1000 times greater than in the fiducial model (set 7, Table [2]). A wide range of cosmic ray 
ionization rates is inferred from models of interstellar chemistry, and observations of HD, OH, D 



medium ( 


Soitzer & Tomasko 


1968: ^ 


kVebber 


1998: 


van Dishoeck &: Black 


1986: 


Le Petit et al. 


2004: 


Federman et al. 


1996: Indriolo et al. 


2007 




McCah et al. 


2003 




Indriolo et al. 


2009). The dis 


tance 


to the ealactic ridge is also important (Melia & Fatuzzo 


2011 


). 



From Fig. [22]we see that ambipolar diffusion loses much of its impact at a cosmic ray ionization 



rate of 10" 



-16 



10- 



-15 



In the most extreme case, the Ohmic and ambipolar dead zones coincide. 



We also varied the dust-to-gas ratio, monomer size and surface density power-law index for a cosmic 
ray ionization rate of 5 x 10"^'' s~^ (not shown here). The ambipolar diffusion still determines the 
dead zone outer edge location and the width of the transitional region. A sharp outer edge was 
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Fig. 20. — MRI-active region appearance for stars of various masses, with luminosities correspond- 
ing to age 1 Myr. Colors and symbols are as in figure SI 
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Fig. 21. — Accretion rate radial profiles in disks around stars with masses from 0.4 to 2 Mq. 
Dotted lines indicate the profiles when the whole disk is active. 
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propose d as a potential pressure trap for solid material in iLvra et al.l (|2009l ): iHasegawa &: Pudritz 
( 2010bl lal). For high cosmic ray ionization rates (two bottom panels in Fig. [22]) . a large region 
with Elsasser number greater than 10 appears. We have a large volume with nearly ideal MHD 
conditions. However, the width of the transitional region is still 1-2 AU, seen from the radial 
separation of the Elsasser number contour lines. The corresponding accretion rates for step- function 
and for smooth function in a(A) are shown in Fig. [23l In the case of smooth a(A), there is a peak 
in accretion rate at 8 AU visible for all models in set 7. At the inner disk boundary, M increases 
from 10"*^ to IQ-'^Moyr-^ 



5. Discussion 

Time Evolution of the Dead Zone: Over the lifetime of a protostellar disk, changes are likely 
in the dust-to-gas ratio fdg as solid material is incorporated into planetesimals and planets. Let us 
consider the possibility that the dust evolves down two paths simultaneously. Most of the aggregates 
remain in place, while a fraction grows in mass through mutual collisions that are "lucky", fast 
enough to overcome electrostatic forces and be compacted, but slow enough to avoid disruption. If 
the monomers' size distribution remains unchanged, the fluffy aggregates' mass fraction relative to 
the gas will decline with time, causing the dead zone to shrink and change shape. In Fig. [2l] left 
panel we show schematically how the dead, transitional and active regions change their midplane 
radial extent with time, and so with decreasing dust-to-gas ratio. Time increases upward and is 
labeled in arbitrary units. The picture is made more complicated by the shorter timescales for 
evolution near the star. 

In our first set of models we studied the impact of the monomer size on the ionization 
state. Real aggregates surely contain monomers wit h a distribution of sizes, possibly resembling 



the interstellar grain size distribution n(a) oc a ( Mathis et al.l 119771 ). Accounti ng for the size 



distribution would likely make the Hall effect stronger at intermediate gas densities (IWardle &: Ng 



19991 ) . Future work should clarify how the dispersion in oq in combination with the radius-dependent 



dust-to-gas ratio changes the dead zone's shape. 

The models in sets 3 and 4 show how changing disk mass and surface density shift the active 
zone inward or outward. In Fig. [24| dashed lines represent the correction to the dead and active 
zone borders when the density slope (or disk mass) changes at the same time the dust-to-gas ratio 
is decreasing. As the outer disk loses mass, the active zone expands inward (compare blue solid and 
dashed lines). The dead zone is situated in the inner part of the disk, where mass is gained slowly 
over time due to the weakness of the MRI-activity in the outer disk. For low dust-to-gas ratio, a 
second MRI-active zone appears inside 10 AU. We expect the accretion flow to drain mass from 
this new zone faster than mass is re-supplied from the outer parts. The surface density may be 
reduced locally, leading to a valley in the surface density. Over the long term the valley in surface 
density near the metal line will grow, becoming better and better ionized and leading to still higher 
accretion rates locally. If this process runs away, it could clear the inner disk of material. The right 
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log(A) tor CR flux 5.d-18 log(A) for CR flux 5,d-17 
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Fig. 22. — Dead, transitional and active zones with the cosmic ray ionization rate varying from 
5 X 10"^® s~^ to 5 X 10^^^ s~^, 5 X 10"^^ s~^ and 5 x 10~^^ s~^. Colors and symbols are as in 
figure m 
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Fig. 23. — Accretion rates vs. radius in set 7 with cosmic ray fluxes at the disk surface from 
5 X 10^^® to 5 X 10~^^. The dotted Une indicates the profile when the whole disk is active. Note 
that the smooth dependence of turbulent stress on Elsasser number leads to higher accretion rates 
in the A > 10 regions (see section 3.2). 
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panel in Fig. [M] shows the expected time evolution of the surface density. Please note that Fig. [23] 
is a schematic representation. A more precise picture requires quantitative disk evolution models. 

Dependence of the Stress on Elsasser Number: Gaps appear in the accretion rate's radial 
profile in many cases when we assume the stress parameter a is a step function across Elsasser 
number unity (figures [TOl [T^ [T8]) . A step function adequately approximates the turbulent stress' 
variation with height in the inner disk, where conductivity gradients are steep. The step function is 
harder to justify when modeling the radial variation in the outer disk, however. We demonstrated 
here that the outer disk can have a quite different accretion rate profile if a (A) is a smooth function. 
The reason is that the Elsasser number is near unity over a large range of radii r > 10 AU. Thus 
much of the outer disk is marginally MRI-active. 

We assume a field strength such that the magnetic pressure is a fixed fraction of the midplane 
gas pressure at each radius. A more complete approac h would inv olve making the stress proportional 



to the magnetic pressure in the saturated turbulence. iBail (120111 ) provides a fitting function for the 
minimum possible plasma beta of the saturated fields in the presence of ambipolar diffusion. We 
have constructed the corresponding turbulent stress profile (not shown in this paper), applying the 
floor for plasma beta mentioned above. It is reassuring that the radial profiles of the accretion rate 
are basically the same as those we present here. When compared to the accretion rates derived 
from a tx 1//3, the accretion is faster near the inner radial boundary using Bai's fit, because the 
weaker magnetic fields leave more of the upper layers with /3 < 1. This deserves to be investigated 
in more detail using 3-D MHD simulations of MRI turbulence including both Ohmic and ambipolar 
diffusion. 



Accretion Variability: It is noteworthy that all models show an accretion rate peak just inside 
the metal line, which lies at 8 AU in our fiducial model. This radial location is affected by the disk 
temperature. The temperature thus affects the location and shape of the MRI-active zone. 

The luminosity produced by the typical T Tauri accretion rate 10~*MQ/yr on a star of Solar 
mass and 2 Solar radii is GMMq/t ~ O.ILq, probably not enough to significantly change temper- 
atures inside the disk region we consider. However during FU Orionis outburs ts, the accretion lumi- 
nosit y rises sufficiently to increase the overall system brightness a hundredfold (jHartmann &: Kenyon 



19961 ). requiring disk temperatures to at least triple. This is likely to expand the Mg"'"-dominated 



region. The outbursts last about 100 years, a few orbits at 10 AU in our fiducial model and long 
enough for MRI turbulence to reach saturation in the newly-coupled gas. Thus we expect the peak 
in the accretion rate to move to larger radii as in set 5. 



Limitations: The outer radial border of the MRI-active zone coincides with the gas density 
contour 10~^^ gcm^'^, lying at 128 AU in our fiducial model (section 4.1). Since gas-phase molecular 
reactions dominate the recombination there and the assumption about fast chemical equilibrium 
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Fig. 24. — Left: Schematic picture of the time evolution of the dead, transitional and active regions. 
Solid lines show the regions' borders when only the dust-to-gas ratio is allowed to decrease with time. 
Dashed lines include corrections for the gas surface density and disk mass evolving simultaneously. 
Right: The black line shows the surface density in our fiducial model. The red line shows the 
expected evolution to form a local gas pressure maximum capable of trapping solid bodies near the 
metal line (8 AU), surrounded by valleys in the surface density which will be better-ionized and 
potentially MRI-active. 
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may become inaccurate, it would be useful to recalculate the ion and electron densities there 
using the extended chemical network, in order to validate our findings on the ambipolar diffusion. 
Additionally, we u se the results from non-stratified simulations of MRI with ambipolar diffusion 
(jBai &: Stonell201lh . to construct our s mooth turbu l ent str ess function and the accretion rates (see 
Section 3.2). In the follow-up work of lBai &: Stond (|2013l ). it is shown that the MRI-active layers 
at lAU can first become turbulent but quickly be stabilized by the ambipolar diffusion. If this 
finding holds, this would affect our predictions of the MRI-driven accretion rate. 

Our chemical network includes only one representative metal, magnesium, but of course many 
metals can occur in the gas phase, and each has its own thermal desorption temperature. There 
should then be multiple metal lines, one for the freeze-out of each metal. 

All the metal lines will shift when we consider a more detailed temperature distribution. In 
our simple model the temperature is constant on cylinders, an approximation which is appropriate 
only for the disk's optically thick inte rior. The surface layers are typically hotter because they are 
directly illuminated by the starlight (jPullemond et al.l 120071 ). so the metal lines would bend away 
from the star with increasing distance from the equatorial plane. 

We may also ask whether a steeper radial variation in the Elsasser number can be obtained 
with a radial gradient in, for example, dust-to-gas ratio or grain charging properties. A steeper 
A(r) would yield bigger differences in the stresses between the transitional, dead and MRI-active 
zones. Such variations also open the possibility of planetesimal traps at multiple radii. 



6. Conclusions 

We determine dead zone sizes and shapes in disks where the ionization state is controlled 
by recombination on fluffy dust aggregates. We study the effects of parameters including the 
constituent monomer size, the aggregates' mass fraction, the magnetic field strength, the masses of 
the disk and star, the slope of the surface density profile, the disk temperature and the cosmic ray 
ionization rate. The three parameters most strongly affecting the dead zone are the monomer size, 
the dust mass fraction and the cosmic ray rate. 

We find that ambipolar diffusion generally defines the dead zone's outer radial boundary, 
as well as the thickness of the transitional layer just inside the boundary. The ambipolar diffu- 
sion is overridden by Ohmic dissipation near the dead zone edge only if the magnetic fields are 
weak enough that the midplane plasma beta is greater than 4000, corresponding to field strengths 
B <0.15 mGauss at 100 AU. 

The fractal aggregates are efficient in sweeping up free electrons and ions, providing a major 
recombination pathway at gas densities greater than 10~^^ gcm~'^. Nevertheless the disk surface 
layers couple to the magnetic fields sufficiently for MRI turbulence to drive accretion. Accurately 
estimating the resulting mass flow rates requires carefully treating the stresses' dependence on the 
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Elsasser number. A step function is an inadequate description because so much of the disk has 
Elsasser numbers near the critical value. The mass flow is W~^-W~^ MQyr~^ at 1 AU for ou r 
fiducial cosmic ray ionization rate of 5 x 10~^^ s~^. This is similar to the findings of iBail ( 201ll ). 
who considered compact grains rather than fractal aggregates. 

The dead zone's outer edge, which we call the transitional region and define by Elsasser numbers 
from 0.1 to 1, is generally very broad at the midplane. Its width is 12 AU in our fiducial model 
and ranges from 2 AU in our most strongly-ionized model from set 7, to the entire disk at the 
maximum dust-to- gas ratio from set 1. The mass flow rates indicate that the transitional region 
will accumulate no sharp density bump as the disk evolves. On the other hand we find that the 
magnetic coupling, and therefore the flow rate, declines steeply where the temperature falls below 
a threshold allowing metal atoms to thermally adsorb on the grains. Such a metal freeze-out can 
cause a local maximum in the flow rate, leading to a gap in the surface density with a long-lived 
pileup of material just outside. This kind of accretion peak occurs in all models across the parameter 
study, shifting in radius according to the disk temperature distribution. 

By varying each of the parameters in turn, we find: 



• At dust aggregate mass fractions near 10"^, thin MRI-active surface layers yield mass fiow 
rates 10~^ MQyr"^. At mass fractions below 10~^ the dead zone takes a fish-tail shape in 
meridional cross-section, and the mass flow rate is around 10~^ MQyr~^. Even at very low 
dust concentrations, there is a patch at r > 10 AU where the magnesium atoms freeze out 
and ambipolar diffusion pushes the Elsasser number into the transitional regime 0.1 < A < 1. 

• The monomer size also affects the dead zone shape. When the monomers are 10 nm or smaller, 
giving the aggregates a large cross-section for a given solid mass, the activity is confined to 
narrow isolated surface layers along with a blob of marginally active gas near 100 AU. When 
the monomers are 100 nm or larger, the dead zone takes the fishtail shape. 

• Varying the surface density slope or the disk mass reveals the dead zone's origins. With 
increasing disk mass we observe larger MRI-active zones in the outer parts, as higher densities 
lead to more frequent collisions with neutrals, thwarting ambipolar diffusion. A similar effect 
comes from reducing the surface density slope so that more of the mass resides far from the 
star. The dead zone's outer edge lies at 20 AU for disk masses between 0.5 and 1.5 times 
the fiducial model. At surface density slopes 1.3 and steeper, ambipolar diffusion isolates a 
midplane marginally-active region outside 20 AU. 

• Making the disk hotter or colder also changes the dead zone's shape. The fishtail appears 
strongest when the magnesium freezeout line lies near the inner edge of the annulus where 
ambipolar diffusion at intermediate heights separates magnetically-active layers in the better- 
ionized surface and denser interior. The fins of the fishtail-shaped transitional zone appear 
to be sensitive to the local ion density and therefor to the disk parameters. 
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• The dead zone's dependence on stellar mass combines the impacts of disk mass, gravity 
and temperature. For our chosen stellar masses and ages (0.4 to 2 Mq, all at 1 Myr), 
the disks are so hot that the Mg+-dominated region extends far from the star. Both the 
magnesium freezeout line and the corresponding accretion rate peak are shifted to larger 
radii with increasing stellar mass. The magnesium line ranges from 10 to 25 AU over the 
stellar mass range we consider. 



We identify a new way to create a local pressure maximum capable of trapping solid particles. 
The steep gradient in mass flow rate across the metal freezeout line means the magnetic stresses 
will evacuate the side of the line near the star, while material piles up beyond. Gas drag forces from 
headwinds outside the maximum, and tail winds inside, will bring solid particles toward the peak. 
The radial gradients could grow strong enough to trigger the formation of vortices, con tributing to 



angu l ar momentuni trari sport in the weakly magnetically coupled transitional region (jFlock et al 



2012; 



Raettig et al.ll2012l ) and further concentrating solid material. Testing the effectiveness of such 



traps is a challenging task. Probing the distribution of the field strength in the MRI turbulence 
requires global 3-D MHD simulations with ambipolar diffusion. Ionization chemistry and the metal 
freezeout must be treated. The pressure profile's evolution near the metal line over the longer 
accretion timescales might be explored using azimuthally-averaged mass transport calculations 
with turbulent stresses obtained from MHD simulations. Since the MRI turbulence and dust are 
mutually coupled and feedback is possible, an important part of the story will be the evolving dust 
abundance and the corresponding changes in the shape of the dead zone. 



N. Dzyurkevich was supported by the Deutsche Forschungsgemeinschaft (DFG) through For- 
schergruppe 759, "The Formation of Planets: The Critical First Growth Phase" . We acknowledge 
support from the Deutsches Zentrum fiir Luft- und Raumfahrt (DLR), support code 50 OR 0401. 
N. J. Turner's participation was supported by the NASA Solar Systems Origins program under 
grant 07-SSO07-0044, and by the Alexander von Humboldt Foundation through a Fellowship for 
Experienced Researchers. The work was carried out in part at the Jet Propulsion Laboratory, 
California Institute of Technology. We thank Prof. W. Brandner and M. Gennaro for the tutoring 
on the STELLAR code. We thank S. Okuzumi for a very thorough referee report that greatly 
improved the quality of this paper. 

REFERENCES 

Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C, & Dullemond, C. P. 2009, ApJ, 700, 1502 

Armitage, P. J., Livio, M., & Pringle, J. E. 2001, MNRAS, 324, 705 

Bai, X.-N. 2011, ApJ, 739, 51 

Bai, X.-N., & Goodman, J. 2009, ApJ, 701, 737 



- 51 - 



Bai, X.-N., & Stone, J. M. 2011, ApJ, 736, 144 
— . 2013, ArXiv e-prints 

Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 
— . 1998, Reviews of Modern Physics, 70, 1 

Beck, R., Brandenburg, A., Moss, D., Shukurov, A., & Sokoloff, D. 1996, ARA&A, 34, 155 
Beck, R., & Hoernes, R 1996, Nature, 379, 47 
Blaes, O. M., & Balbus, S. A. 1994, ApJ, 421, 163 

Blum, J. 2004, in Astronomical Society of the Pacific Conference Series, Vol. 309, Astrophysics of 
Dust, ed. A. N. Witt, G. C. Clayton, & B. T. Draine, 369 

Blum, J., & Wurm, G. 2000, Icarus, 143, 138 

Blum, J., Wurm, G., Poppe, T., & Heim, L.-O. 1998, Earth Moon and Planets, 80, 285 

Bodenheimer, P., Laughlin, G. P., Rozyczka, M., & Yorke, H. W., eds. 2007, Numerical Methods 
in Astrophysics: An Introduction 

Bouwman, J., Meeus, G., de Koter, A., et al. 2001, A&A, 375, 950 

Brandenburg, A., Nordlund, A., Stein, R. F., k Torkelsson, U. 1995, ApJ, 446, 741 

Brauer, F., Henning, T., & Dullemond, C. P. 2008, A&A, 487, LI 

Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368 

Cowling, T. G. 1976, Magnetohydrodynamics (Monographs on Astronomical Subjects, Bristol: 
Adam Hilger, 1976) 

Desch, S. J. 2004, ApJ, 608, 509 

Dominik, C, Blum, J., Cuzzi, J. N., & Wurm, G. 2007, Protostars and Planets V, 783 
Draine, B. T., Roberge, W. G., & Dalgarno, A. 1983, ApJ, 264, 485 

Dullemond, C. P., Hollenbach, D., Kamp, I., Sz D'Alessio, P. 2007, Protostars and Planets V, 555 
Dzyurkevich, N., Flock, M., Turner, N. J., Klahr, H., & Henning, T. 2010, A&A, 515, A70 
Federman, S. R., Weber, J., & Lambert, D. L. 1996, ApJ, 463, 181 
Flaig, M., Ruoff, P., Kley, W., & Kissmann, R. 2012, MNRAS, 420, 2419 
Fleming, T., & Stone, J. M. 2003, ApJ, 585, 908 



- 52 - 



Flock, M., Dzyurkevich, N., Klahr, H., Turner, N., & Henning, T. 2012, ApJ, 744, 144 

Flock, M., Dzyurkevich, N., Klahr, H., Turner, N. J., & Henning, T. 2011, ApJ, 735, 122 

Fromang, S., & Nelson, R. P. 2006, A&A, 457, 343 

Furlan, E., Watson, D. M., McClure, M. K., et al. 2009, ApJ, 703, 1964 

Gammie, C. F. 1996, ApJ, 457, 355 

Glassgold, A. E., Najita, J., & Igea, J. 1997, ApJ, 480, 344 
Guan, X., & Gammie, C. F. 2011, ApJ, 728, 130 

Hartmann, L., Calvet, N., Gullbring, E., & D'Alessio, P. 1998, ApJ, 495, 385 
Hartmann, L., & Kenyon, S. J. 1996, ARA&A, 34, 207 
Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJS, 82, 167 
Hasegawa, Y., & Pudritz, R. E. 2010a, ApJ, 710, L167 
— . 2010b, MNRAS, 401, 143 
— . 2011, MNRAS, 417, 1236 

Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742 
Hawley, J. F., & Stone, J. M. 1998, ApJ, 501, 758 
Hirose, S., & Turner, N. J. 2011, ApJ, 732, L30 

Hughes, A. M., Wilner, D. J., Qi, C., & Hogerheijde, M. R. 2008, ApJ, 678, 1119 

Igea, J., & Glassgold, A. E. 1999, ApJ, 518, 848 

Ilgner, M., & Nelson, R. P. 2006, A&A, 445, 223 

Indriolo, N., Fields, B. D., & McCall, B. J. 2009, ApJ, 694, 257 

Indriolo, N., Geballe, T. R., Oka, T., & McCall, B. J. 2007, ApJ, 671, 1736 

Johansen, A., Klahr, H., & Henning, T. 2011, A&A, 529, A62 

Johansen, A., & Youdin, A. 2007, ApJ, 662, 627 

Johansen, A., Youdin, A., k Klahr, H. 2009, ApJ, 697, 1269 

Kempf, S., Pfalzner, S., & Henning, T. K. 1999, Icarus, 141, 388 

Kim, W.-T., & Ostriker, E. C. 2000, ApJ, 540, 372 



- 53 - 



Krejny, M., Matthews, T. G., Novak, G., et al. 2009, ApJ, 705, 717 
Kretke, K. A., & Lin, D. N. C. 2007, ApJ, 664, L55 
Kunz, M. W., & Balbus, S. A. 2004, MNRAS, 348, 355 
Le Petit, F., RouefF, E., & Herbst, E. 2004, A&A, 417, 993 

Lyra, W., Johansen, A., Zsom, A., Klahr, H., & Piskunov, N. 2009, A&A, 497, 869 
Mac Low, M.-M., Norman, M. L., Konigl, A., & Wardle, M. 1995, ApJ, 442, 726 
Martin, R. G., k Livio, M. 2012, MNRAS, 425, L6 

Martin, R. G., Lubow, S. H., Livio, M., & Pringle, J. E. 2012, MNRAS, 423, 2718 

Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 

McCall, B. J., Huneycutt, A. J., Saykally, R. J., et al. 2003a, Nature, 422, 500 

— . 2003b, Nature, 422, 500 

Meakin, P. 1991, Reviews of Geophysics, 29, 317 

Meha, F., & Fatuzzo, M. 2011, MNRAS, 410, L23 

Miller, K. A., & Stone, J. M. 2000, ApJ, 534, 398 

Muranushi, T., Okuzumi, S., &: Inutsuka, S.-i. 2012, ApJ, 760, 56 

Okuzumi, S. 2009, ApJ, 698, 1122 

Okuzumi, S., & Hirose, S. 2011, ApJ, 742, 65 

Perez-Becker, D., & Chiang, E. 2011, ApJ, 727, 2 

Pinilla, P., Birnstiel, T., Ricci, L., et al. 2012, A&A, 538, A114 

Raettig, N., Lyra, W., & Klahr, H. 2012, ArXiv e-prints 

Sano, T., Inutsuka, S.-i., Turner, N. J., & Stone, J. M. 2004, ApJ, 605, 321 

Sano, T., & Miyama, S. M. 1999, ApJ, 515, 776 

Sano, T., Miyama, S. M., Umebayashi, T., & Nakano, T. 2000, ApJ, 543, 486 
Sano, T., & Stone, J. M. 2002a, ApJ, 570, 314 
— . 2002b, ApJ, 577, 534 

Semenov, D., Wiebe, D., & Henning, T. 2004, A&A, 417, 93 



- 54 - 



Simon, J. B., Bai, X.-N., Stone, J. M., Armitage, P. J., & Beckwith, K. 2012, ArXiv e-prints 

Simon, J. B., Hawley, J. F., & Beckwith, K. 2009, ApJ, 690, 974 

Spitzer, Jr., L., & Tomasko, M. G. 1968, ApJ, 152, 971 

Suyama, T., Wada, K., & Tanaka, H. 2008, ApJ, 684, 1310 

Tout, C. A., & Pringle, J. E. 1992, MNRAS, 259, 604 

Turner, N. J., Carballido, A., & Sano, T. 2010, ApJ, 708, 188 

Turner, N. J., & Drake, J. F. 2009, ApJ, 703, 2152 

Turner, N. J., & Sano, T. 2008, ApJ, 679, L131 

Umebayashi, T. 1983, Progress of Theoretical Physics, 69, 480 

Umebayashi, T., & Nakano, T. 1981, PAS J, 33, 617 

— . 2009, ApJ, 690, 69 

Uribe, A. L., Klahr, H., Flock, M., & Henning, T. 2011, ApJ, 736, 85 

van Dishoeck, E. F., & Black, J. H. 1986, ApJS, 62, 109 

Wardle, M. 2007, Ap&SS, 311, 35 

Wardle, M., & Ng, C. 1999, MNRAS, 303, 239 

Wardle, M., & Salmeron, R. 2012, MNRAS, 422, 2737 

Webber, W. R. 1998, ApJ, 506, 329 

Weidenschilling, S. J., & Cuzzi, J. N. 1993, in Protostars and Planets III, ed. E. H. Levy &: J. I. 
Lunine, 1031-1060 

Wurm, G., & Blum, J. 1998, Icarus, 132, 125 

Zhu, Z., Hartmann, L., k Gammie, C. 2010, ApJ, 713, 1143 



This preprint was prepared with the AAS lATp^X macros v5.2. 



- 55 - 



A. Determining Gas-Phase Metal Abundance 

Most of the metal atoms are incorporated into grains. Tiie remainder can occur either free 
in the gas, where they are available to participate in the ionization balance, or adsorbed on grain 
surfaces. We determine here the number density of the metal atoms in the gas phase, which depends 
on temperature. We consider only the atoms' thermal adsorption on and desorption from the grains. 
Gas-phase species X when captured on the grains is labeled gX. The equilibrium number density 
on the grains is given by 

n[gX] = -l^n[X]ng, = an[X], (Al) 

where the coefficients ki = Ugrfth and k2 are adsorption and desorption rates. The number density 
in the gas 

n[X] = n[X]*°'^l - n[gX]. (A2) 
We determine the amount of magnesium on the grain surfaces using the expression 

„bMg| = „:i^, (A3) 

1 + a 

and in the gas phase 

n[Mg] = . (A4) 

1 + a 

The initial ab undances are 1 for Ho , 9.75x10"^ for He, 3.62 and 8.53 xlO""^ for C and O, and 1 x 
10~^^ for Mg ( Ilgner fc Neisonlbnnd v Not every collision with a grain leads to sticking, and the stick- 



ing probability is Sx = exp{-Ekin/ \J {2ED^Etrans)), with AEtrans '^^^^ 2 x 10 ^ eV for the energy 
transfer to the grain via lattice vibrations. E^^ = 450, = 1510 and E^^ = B SOOkn K. The 



desor ption coefficient /c2 = fo ^^pi—Eo/kT) with i^q = a/(3 x lO^^cm "^E^) / {ir'^'mx)) ( Hasegawa et al 



19921 ). After finding the gas-phase metal abundance, we compute the total recombination rate q 



for the dominant ion by applying an iterative approach (section [2 

The transition between two representative ions leads to the jump in Elsasser number and in 
the corresponding accretion rates (Fig. [25]) . Here we show how the accretion rates would look like, 
if the representative ion is chosen to be Mg+ everywhere, or HCO"*" everywhere. The accretion 
rates for the fiducial model follow the case of pure Mg"*" inside of r < 5AU, and the case of pure 
HCO+ outside of r > 9AU. 



B. Gaussian Charge State Distribution — Limit of Validity 



We find the abundances of the ions, electrons and charged grains from lOkuzumil (j2009l ) eqs. 32, 
33 and 22, letting the grains charge by up to 30 electrons either side of their mean charge. The 
grains' charge state distribution is very well-approximated by the Gaussian when the grains carry 
many charges each. Small aggregates however charge weakly enough that their few discrete charge 
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states can differ noticeably from a Gaussian. As a result the plasma is not quite charge-neutral 
overall. Fig. [26] demonstrates the accuracy of the charge neutrality for three different radii in 
the fiducial model. The largest errors occur near the midplane where most of the charge resides 
on the grains. Accuracy is good when the fluffy grains are made of A'^ > 400 monomers, while 
relatively large departures from neutrality occur when N < 10. Several of the assumptions used by 



Okuzumil (j2009l ) break down when the number of monomers is too small, including the assumption 



of negligible electric polarization. But most importantly, charge neutrality is violated when the 
grains' charge spread is too small to be represented with a Gaussian. 



C. MRI Criterion When Ambipolar Diffusion Dominates 

Whether MRI turbulence can operate under an ambipolar-dominated magnetic diffusivity is 
sometimes judged using the ion magnetic Reynolds number. In this Appendix we examine the 
general conditions for this number to be equivalent to the ambipolar Elsasser number. That is, we 
seek the requirements for Rci A^ or 

or 

riA^j • (C2) 



fap 1 



Recall 

from eq. [27l Let each of the three conductivities be dominated by a single species. Label the 
main Ohmic species' mass niQ, charge qo, number density no, collision coefficient 70 and coupling 
parameter bo, and similarly for the Hall and Pedersen conductivities: 

ec f eBqo \ , . 

<yo = -^noqo (C4) 

ec nnqn 
^ 1 + ( ^Bqn 



= -5 ; — I — 7T y^^) 



B 1 

' Xmpc^pp 



We have used charge neutrality to derive eq. IC5I from eq. [241 Notice that ao is independent of 
the field strength, so we can end up with rjA oc B^ as in eq. IC2I only if the final term in eq. IC3I is 
negligible. 



Now the perpendicular conductivity a± = y crj^ + ap can be dominated by either an or ap. 
We consider the two possibilities in turn. If an ^ crp, the remaining piece of the right-hand side 



-57- 



step function «(A) 




smooth function a(A) 

IQ-^F . . . ■ ■ ■ r- 




Fig. 25. — Accretion rates in solar masses per year for models with the representative ion Mg"*" or 
HCO"'' (red and purple lines). Black line shows the fiducial model with the transition between the 
two ions, which occures between 5 and 9 AU. The dotted line shows the accretion rate for the case 
when the whole disk is active with the turbulent stress from eq. 29. Note that the smooth-function 
turbulent stress can lead to higher accretion rates in the A > 10 regions (see section 3.2). 




Fig. 26. — Left: Charge neutrality for R = 2, 5 and 9.4 AU, for N = 400. Right: same labeling for 
charge neutrality at the midplane, for variable number of monomers N. 



- 58 - 



of eq. 



IS 



'H 



eBqp 
mpc^pp 



This reduces to eq. IC2l if three conditions are met: 



1 + 



(C7) 



1. The ions dominate the Pedersen conductivity. 

2. The ions have the same number density and charge magnitude as the Hall species, i.e. 
hh = Up = Tii and {qnl = \qp\ = \qi\- For example, the ions could also dominate the 
Hall conductivity. 



3. Both the ions and the Hall species are poorly coupled, so bi 
terms not quadratic in B can be neglected. 



bp <^ 1 and 6/^ <C 1 and the 



In the other case, where cth <^ crp, the right-hand side of eq. IC3I becomes 



1 

ap 



npqp 




mpjpp 



eBqp 
mpc^pp 



(C8) 



This reduces to eq. IC2I if the ions again determine the Pedersen conductivity, and in addition the 
term proportional to B'^ dominates, i.e. the ions couple to the fields and bi = bp ^ 1. 

We next look at how the Ohmic term becomes negligible in the ambipolar diffusivity when the 
other requirements for eq. IC2I are met. In the first case, where an ^ o"p, the Ohmic term is small 
if ao ^ ^h/'^p or 



noqobo > 



l + bl 



(C9) 



Now in this first case, eq. 
eq. IC9I simplifies to 



\npqpj bpil + bjjf 
works only when conditions 1-3 are satisfied. Under the conditions. 



bo > 



riiC 



(CIO) 



noqo 

That is, the Ohmic species must not decouple too much from the magnetic fields. Throughout our 
fiducial model disk, and over a wide range of other models, the Ohmic conductivity comes from 
the electrons. Since these have the same charge magnitude as the ions, and a coupling parameter 
bf. ~ lOOOftj due to their low mass, eq. ICIOI means the ion coupling must satisfy 10006? ^ iii/rif.. 
Also the electrons have the same number density as the ions in the ion-electron plasma limit where 
grain-surface recombination is unimportant, while in the opposit e limit of an ion -dust plasma the 



electrons' depletion onto grains yields rie « lO-^n^ (eq. 38 of lOkuzumil l2009l ). Eq. [CTO] then 



requires bi S> 0.03 and bi 3> 0.3, respectively. These are incompatible or just barely compatible 
with bi <^ 1 needed by condition 3. Eq. IC2I is therefore valid at most in a narrow range of ion 
coupling parameters when an ^ ap. 
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In the second case, where an ^ cp, the Ohmic term can be dropped when ao ^ crp, that is 
when 

npqpbp 



noqo^O > . ,2 ■ (Cll) 

i \ Op 

Since eq. IC2I works in this case only if fep ^ 1, it is enough to say the Ohmic term is smaU when 



ho » ^hT\ (C12) 
noqo 

the same form as eq. IClOi Even the stronger of the two constraints derived above for eq. IClOt 
hi ^ 0.3, is aheady met in this case by our requirement that 6j = 6p ^ 1. 

Considering the two cases together, we see that the ion Reynolds number is equivalent to the 
ambipolar Elsasser number if the ions couple to the field and dominate the Pedersen conductivity, 
which in turn is larger than the Hall conductivity. 
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